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We describe the oligosaccharides-exchange dynamics performed by so-called D-enzymes on polysaccharides. 
To mimic physiological conditions, we treat this process as an open chemical network by assuming some of 
the polymer concentrations fixed (chemostatting). We show that three different long-time behaviors may 
ensue: equilibrium states, nonequilibrium steady states, and continuous growth states. We dynamically and 
thermodynamically characterize these states and emphasize the crucial role of conservation laws in identifying 
the chemostatting conditions inducing them. 


PACS numbers: 05.70.Ln, 05.70.-a, 82.20.-w 

I. INTRODUCTION 

Biological systems use large and branched chains of 
basic sugars, called polysaccharides, to store energy^. 
Glucans such as glycogen and starch are polysaccharides 
whose building blocks are D-glucose monosaccharides. 
Despite the apparent simplicity of their constituents, 
their metabolism involves several chemical steps, each 
performed by a specific set of enzymes^. Interestingly, 
some of these catalysts lack specificity regarding the re¬ 
action they catalyze or the substrates they act on^. An 
example is provided by (1^4)-alpha-D-glucans^’^^^ (EC 
2.4.1.25), also called D-enzymes, which act on pairs of 
glucans regardless of their size ' . Specifically, D-enzymes 
catalyze the transfer of groups of glycosyl residues from 
a donor glucan to an acceptor glucan^’®. Experimental 
evidences highlight the presence of bonds between gly¬ 
cosyl residues which are not cleaved by D-enzymes^—at 
least not over physiological time scales ' . These bonds 
are called forbidden linkages^. In this way, D-enzymes 
transfer segments of glucan chains containing one or more 
forbidden linkages, and the transfer of segments contain¬ 
ing one forbidden linkage are the most probable^. Also, 
each glucan chain is characterized by a reducing-end glu¬ 
cose which is not transferred by D-enzymes"^’^. Hence, 
glucans made of just the reducing end can act only as 
acceptor in the transfer. 

Qualitatively, D-enzymes process medium-size glucans 
by disproportionating them into unit-size and big-size 
glucans®. Since their transfers reactions are neutral 
energetically^’®, entropy is the main driving force in this 
system. In closed conditions, this system evolves towards 
an equilibrium state maximizing the entropy^’®. 

In this paper we consider a simplified kinetic description 
of the D-enzyme’s action on glucans, which we treat as a 
chemical network. Since metabolic processes should be 
thought of as part of an open system continuously fed from 
the environment, we mimic these physiological conditions 
by introducing chemostats (i.e. species whose concentra¬ 
tions are kept constant by the environment). Our goal 
is to characterize the dynamical and thermodynamical 


implications of treating the action of the D-enzymes on 
glucans as an open chemical network. In the framework 
of deterministic chemical networks endowed with mass 
action kinetics, we prove that chemostatting can induce 
three different types of long-time behaviors: equilibrium, 
non-equilibrium steady state, and continuous growth. The 
equilibrium state corresponds to the stationary concen¬ 
tration distribution in which the concentration currents 
along each reaction pathway vanishes (detailed balance 
property^®). Non-equilibrium steady states refer to sta¬ 
tionary distributions violating detailed balance. Hence, 
contrary to equilibrium states, a continuous and steady 
flow of mass circulates across the network. Finally, the 
continuous growth regime we observed corresponds to 
a non-stationary state characterized by continuous and 
steady flow of mass entering the network, and resulting in 
its continuous growth. We emphasize the dynamical and 
thermodynamical role of conservation laws and emergent 
cycles in identifying the chemostatting conditions leading 
to these states. We are thus able to confirm the general 
relation between the number of chemostatted species and 
the number of independent thermodynamical forces—or 
affinities —found in Ref.^^. Despite the simplicity of our 
description, the closed system results found in Ref.” are 
reproduced and the qualitative disproportionating behav¬ 
ior of D-enzymes® is captured by our (chemostatted) open 
system description. 

The plan of the paper is as follows: in sec. H the kinetic 
model is established and the related rate equation descrip¬ 
tion for the concentration of polysaccharides is introduced. 
In sec. HI the chemostatting conditions leading to non¬ 
equilibrium steady states rather than equilibrium ones 
are found. For this purpose, both the conservation laws of 
the dynamics and the emergent cycles of the network are 
analyzed. The dissipation of the non-equilibrium steady 
state is also studied. The network’s conservation laws 
identified in sec. HI A are used in sec. IV to derive the 
steady-state concentration distributions for different num¬ 
bers of chemostats. The explosive asymptotic behavior 
is described in sec. V. Conclusions are drawn in sec. VI. 
Some technical derivations and proofs are provided in the 
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FIG. 1: (a) The typical monomer-exchange reaction de¬ 
scribing the action of D-enzymes on glucan chains, (b) 
The attachment of free monomers to other species is not 
allowed. 


appendices. 


II. THE KINETIC MODEL 

The action of D-enzymes is modeled as follows (see also 
Fig. 1). Glucans are treated as polymers whose monomers 
represent single transferable segments. Hence, each glucan 
is identified by its number of monomers, or equivalently 
by its monomeric mass. The enzymatic steps performed 
by the D-enzymes in order to achieve the transfer are not 
explicitly described—they are coarse-grained—, and we 
describe the interaction between two polymers of mass n 
and m as a mass-exchange process^^: 

(n) -I- (m) {n1) {m - 1), 

for n > 1, m > 2, (1) 

where Knm denotes the related coarse-grained rate con¬ 
stant. Transfers of oligosaccharides longer than one 
monomeric unit are less probable^, and are not considered 
in our description. We take into account the presence of 
non-transferable units by imposing the size of the donor 
glucan (m) to be greater than one^’^. 

Let us note that each reaction is reversible because the 
backward path is already included in (1) (it is realized 
by replacing n —>■ m — 1 and m —> n -I- 1 in the above 
expression). Furthermore, the constraint on the minimal 
size of the donor molecules imposes that m > 2. Since 
we describe the glucans as linear polymers, and since 
D-enzymes do not discriminate the size of the polymers, 
we assume a constant kernel for the reactions: Kmn = 

Vn > 1, Vto > 2. This assumption is based on the evidence 
that the free-enthalpy release resulting from any reaction 
is almost vanishing^’®. Indeed, for any bond cleaved, a 
new one of the same kind will be formed. 

Assuming a large and well stirred pool of interacting 
polymers, the evolution of the system is well described by 
reaction rate equations^^. According to this mean-field 
description, the molar concentration of polymers of mass 
k at time t, = Z^(t), satisfies the following first order 
differential equations 

i E <rn , for A: > 1. (2) 

n>l ' --^ 

m>2 ^J<^m 


The ^ factor in front of the summation takes into account 
that summing over all n > 1 and m > 2 includes every 
reaction pathway twice^®. represents the element of 

the stoichiometric matrix related to the species of mass 
k and to the reaction involving an acceptor and a donor 
polymer of mass n and m, respectively. The reaction 
scheme in (1) implies that 

VL = t+l + ^1-1 - - <^1, (3) 

where 5j represents the Kronecker delta. Assuming a 
mass action kinetics, the forward (denoted by -I-) and the 
backward fluxes (—) can be written as: 

7+”"* = , J-""* =(4) 


where Z'^, denotes the concentrations of the polymers 
of size n. To simplify the following discussion, we will 
use the Einstein summation notation: upper indexed 
quantities represent vectors, lower indexed ones covectors, 
and repeated indexes implies the summation over all the 
allowed values for those indexes (1 < n < nmax and 
2 < m < rrimax, or 1 < fc < fc^ax, where rimax, Wmax 
and fcmax are finite in closed systems but infinite in open 
ones). To avoid confusion, exponents will always act on 
parenthesis (e.g. (a)" denotes the quantity a to the power 
n). 

The rate equations (2) assume the following form when 
the expressions for both the stoichiometric matrix (3) and 
the fluxes (4) are considered 

= kZ {Z^ - Z^) -H kZ^Z\ 

Z^ = kZ (Z'=+i - 2Z'= -h -h (5) 

-H kZ^ - Z'=-i) , for fc > 2, 


where Z = denotes the total concentration. 

The second term in the right hand side of (5) arises 
from the constraint that the donor species cannot be 
monomers^^ (see Fig. lb). 

To model the open system we now assume that the 
environment keeps the concentrations of some species 
constant by refilling the consumed ones and eliminating 
the produced ones, see Fig. 2. We call these species 
chemostats^^ and we denote them with the indices ky G 
Hyj where Hy C N represents a subset of all species. The 
remaining (variable) species are explicitly denoted by k^- 

By definition, the chemostats’ concentrations must 
remain constant, Z^y = 0. The rate of chemostatted 
molecules consumed by the reactions in the network 
must therefore be balanced by the rate of chemostatted 
molecules injected/rejected from the system. The rate of 
injection/rejection of the fcy-th chemostat is quantified by 
the external currents^^, whose expression is 


jky _ ^j+nm _ j—nm'^ _ 

= kZ {Z"^ - Z^)-G kZ^Z^ ifky = l 

= kZ ( z '=^+^ - 2Z^y + 

+ kZ^ {Z'^y - Z’^y-^) if ky > 2. 
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FIG. 2: Pictorial representation of a reaction involving a 
chemostat. When a reaction produces a chemostat (here 
a dimer), the environment extracts one molecule of this 
species from the system (dotted light green reaction). On 
the other hand, when a chemostat reacts, a new molecule 
is injected into the system (dashed dark green reaction). 


III. STEADY STATES: CONSERVATION LAWS, 
CYCLES AND DISSIPATION 


Three different types of long-time behaviors have 
been identified for our kinetic model: equilibrium, non¬ 
equilibrium steady state and continuous growth. We start 
by focusing on the chemostatting conditions leading to 
equilibrium or non-equilibrium steady states. The exis¬ 
tence and uniqueness of the steady state is currently a 
priori assumed. 

Closed systems always reach an equilibrium steady 
state^® defined by = 0,Vfcx and = 0,Vn,TO. 
Their dynamics is constrained by conservation laws^^’^'^’^^, 
which fully characterize the equilibrium concentration dis¬ 
tribution. Chemostatting generic chemical species may 
break these conservation laws and may create chemical 
forces—also called ajfinities^^. The appearance of affini¬ 
ties is directly related to that of so-called emergent cycles, 
through which the external chemical forces can act. In 
finite chemical networks, if no emergent affinity arises 
from the chemostatting procedure, the system will al¬ 
ways relax to a unique equilibrium state compatible with 
the chemostats and the non-broken conservation laws^^’^®. 
When emergent cycles—or equivalently affinities—are gen¬ 
erated, the system may evolve towards a non-equilibrium 
steady state defined by = 0, Vfcx and J""* ^ 0 (non¬ 
equilibrium steady state quantities are denoted by an 
overbar in the text). In the following subsections we 
analyze how the closed system’s conservation laws and 
emergent cycles are modified by the gradual increase of 
the number of chemostatted chemical species. In the last 
subsection we relate these to the dissipation in the system. 


A. Conservation Laws 


Conservation laws denote the presence of physical quan¬ 
tities which are conserved during the evolution of the 
system, the so-called components. In general, they can 
be identified from the cokernel space of the stoichiomet¬ 
ric matrix^^d7,i8^ Indeed, if Ik G coker V, namely if 


Ik^nm = 0, the scalar IkZ^ is conserved 


- {ikZ^) = IkZ^ 


= 2 - J-"^) = 0 . 


(7) 


For the closed system, the equation leading to the 
conservation laws is l^+i ~ ~ ^ ^ n < 

^max “ ^max 1 and 2 ^ 7Ti ^ rn^na-K. — ^max ■ If exhibits 
the following solutions: = a and = a ■ k (where 

a is an arbitrary constant, which is taken as one when 
expressing the components), which correspond to the 
conservations of the total concentration Z = Z 

and the total mass M = J2k=T kZ^, respectively. Hence, 
fcmax = M — Z 1. 

However, when the system is opened by setting 
chemostats, the relevant stoichiometric matrix becomes 
the stoichiometric submatrix of the variable species: 

Also, fcmax = oo. No matter what the sizes of the 
chemostatted glucans are, neither the total concentration 
conservation law = a nor the total mass conservation 
law = afcx survives (i.e. they are not anymore ele¬ 
ments of the cokernel space of V^^). We therefore say 
that the total mass and the total concentration are broken 
conservation laws. Nevertheless, when just one chemostat 
is present, fly = {^y}i a new conservation law emerges: 


4!^ = a (fcx - ky). 


( 8 ) 


Hence, the system exhibits just one (net) broken conser¬ 
vation law. It corresponds to the component 


q = M- kyZ, (9) 

which can assume any value in M and takes into account 
that the total mass can change in the system only by 
multiples of the chemostat mass, ky. In presence of more 
than one chemostat, no conservation law survives. 

The components derived in this section —M and Z 
for the closed system and q for the network with one 
chemostat—will be used to characterize the steady state 
distribution in sec. IV. 


B. Emergent cycles 

A cycle represents a finite set of reactions which leave 
the state of the network unchanged. Algebraically they are 
represented as vectors c""* and they belong by definition 
to the kernel space of the stoichiometric matrix (c”"* G 
kerV): iV^„c"- = 0. 

The steady-state currents satisfy J""* = 0 and can 
always be written as linear combinations of cycles. The 
cycle space of our polymers system is however infinite di¬ 
mensional and its complete characterization is of little use. 
However, in order to characterize non-equilibrium steady 
states only the emergent cycles —those cycles that may 
appear when chemostatted species are introduced—are 
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needed^^. Physically, they represent cyclic transforma¬ 
tions leaving the variable species fcx unchanged, but which 
would change the concentrations of the chemostats ky if 
they were not kept constant and contribute to the external 
currents. 

An emergent cycle (y"™) is thus defined by 
= 0. 

1 (10) 
^ 0 for at least one ky, 


C. External Currents and Dissipation 

We now show that at steady state, the emergent cy¬ 
cles determine the external currents and the entropy 
production rate S. 

We first observe that the steady-state external currents 
are in general linear combination of the coefficients 
iyj[ and must satisfy the same constraints (eq. 12). Indeed, 
the steady-state equations in presence of chemostats 


where < !■ denotes the amount of chemostats 

V J ky^Cl'Y 

of mass ky injected (minus sign) or rejected (plus sign) 
from the chemical network during the transformation y"™. 
These quantities cannot take arbitrary values, due to the 
constraints imposed by the conservation laws of 
Indeed, for any conservation law, Z), , a constraint of the 
following form holds 
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2 nm I 


= 0 . 


( 11 ) 


Taking into account the total concentration = a 
and total mass Z). ^ = ak conservation laws, derived in 
sec. Ill A (the emergent conservation law Z), Ms a linear 
combination of the first two on the whole set of species 
indexes), we obtain the following constraints 


Sfcy kyUj^ = 0. 


( 12 ) 


Non-trivial solutions of this set of equations signal the 
presence of emergent cycles, and thus of independent 
affinities, which read'^^ 

( 13 ) 

nm ky 


The set of linearly independent solutions of (12) gives the 
number of independent emergent cycles in the chemostat- 
ted chemical network. If we normalize this set so to have 
the smallest non-vanishing integer values for , these 
values indicate the number of chemostatted species which 
are introduced in or rejected from the system in precisely 
one (emergent) cyclic transformation. 

For less than three chemostats, only trivial solutions of 
(12) exist and therefore no emergent cycle appears. For 
three chemostats, we obtain one emergent cycle charac¬ 
terized by the following normalized values for : 

= kys - ky2, 

= ftyl - fcy3, (14) 

= ky2 - kyi , 

where fcyi, ky2 and ky^ represent the masses of the 
chemostats. For any additional chemostat we obtain 
an additional emergent cycle, each characterized by its 
value for the coefficients v^y. 


iyfex 
2 nm 


-V 

2 ^ nm 


jnm 
ky jnm 



(15) 


are equivalent to eq. (10): the emergent cycles y”"* are 
substituted by the steady state currents J”™ and the 
coefficients v^y by the steady-state external currents I^y. 
Thereby, if no cycle emerges d_ue to the chemostats, the 
steady-state external currents I^y are vanishing, provided 
that the steady state exists. The system is then at equi¬ 
librium. 

The dissipation at steady state is intimately related 
to the external currents^^. Indeed, the (non-negative) 
entropy production rate for our chemical reaction network 
can be written as 


1 T+nm 

S = -V J”™i?ln-- = 

2 / J—nm 

nm 

= -Vz'=^i?ln^-V/'=^i?ln^, ( 16 ) 




= Sx 


where R is the gas constant. At the steady state, the 
internal species’ contribution Sx always vanishes. Hence, 
the dissipation is characterized by the contribution due 
to the chemostats Sy, which is non-vanishing if the set of 
steady-state external currents I^y is also non-vanishing. 
We also mention that, the steady state entropy production 
can be expressed as the sum along a set of independent 
emergent cycles of products of affinities (13) and emergent 
cycle currents^^ Jy-. E = 

Summarizing, the conservation laws provide us with 
both the components— which are useful for expressing the 
steady state distributions—and the constraints (eq. 12) on 
the emergent cycles of the network (eq. 10). Due to these 
constraints, the first emergent cycle appears in the system 
with three chemostats. For any additional chemostat an 
additional independent cycle emerges. Through these 
cycles the environment exerts chemical forces, which are 
generated by the chemostats concentrations. The above- 
defined external currents result from these forces and 
characterize the dissipation. 

We emphasize that the relation between the number 
of chemostats s^, of net broken conservation laws b, and 
of emergent cycles a, is in perfect agreement with the 
general result obtained for finite-dimensional phase space 
in Ref.^^ stating that 

= 6 -I- a 


(17) 
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number of broken independent asymptotic 

chemostats, c. laws, b affinities, a behavior 


0 

0 

0 

ES 

1 

1 

0 

ES 

2 

2 

0 

ES/growth 

3 

2 

1 

NESS/growth 

4 

2 

2 

NESS/growth 


TABLE I: Summary of the behaviors of our model for dif¬ 
ferent numbers of chemostats (ES stands for “equilib¬ 
rium state” whereas NESS for “non-equilibrium steady 
state”). The number of broken conservation laws and in¬ 
dependent affinities are also reported. The growth state 
occurs whenever the concentration of the largest chemo- 
stat is larger than the concentration of the smallest one: 

larger > 


These results are summarized in Tab. I. 



polymer mass, k 


FIG. 3: Equilibrium concentration distribution for the 
closed system of monomers-exchanging polymers at dif¬ 
ferent values of the total concentration Z and total mass 
M. The dark blue bar plot refers to the choice Z = 10 and 
M = 15, while the light blue one to Z = 10 and M — 55. 


IV. THE STATIONARY DISTRIBUTIONS 


We now use the components introduced in the previous 
section to derive the steady-state concentration distribu¬ 
tion for different number of chemostats. The conditions on 
the chemostats’ concentrations not leading to the steady 
state solution are also identified. 

From the steady-state equations corresponding to (5) 
and from the equations for the external currents (6), we 
can write a general expression for the steady-state concen¬ 
trations as a function of the concentration of monomers, 
Z^, the fraction of polymers larger than monomers, 
f = 1 — Z^/Z, and the chemostats fluxes, 




+ E 

ky G Y 


/fcy 

K 


1 — {rY 

1 — r 


Q {k — ky 


1 ), ( 18 ) 


where 0(-) represents the discrete step function (we refer 
the reader to appendix A for details). Here, the number of 
chemostats is arbitrary, and since the external currents at 
steady state satisfy the same constraints as in (12), only 
— 2 of them are independent. In the next paragraphs 
we will discuss in detail the above expression for zero, one, 
two and three chemostats, and the variables Z^, r and 
I^y will be expressed in terms of the components and of 
the chemostats’ concentrations. 


A. Closed system 

As previously discussed, the closed system exhibits 
the following components: Z = J2k=T ^ — 

J2k=T kZ’^. In order to express the equilibrium distribu¬ 
tion algebraically as function of Z and M we consider 
the following limit M ^ Z. In this way fcmax ~ oo and 


imposing Z = ^ 

pression (18) we can write Z^ and r as functions of Z 
and M. Hence 


Z^ = 


M \ MJ 


(19) 


Fig. 3 shows the behavior of this distribution for different 
values of Z and M. As expected, the higher the ratio 
between the mass and the concentration M ^ Z, the 
broader the distribution. 

Remark. The equilibrium distribution we obtained 
from our dynamical description is equivalent to the result 
obtained using maximum entropy approaches and is con¬ 
sistent with experimental observations^. The equivalence 
is inferred by comparing eq. (19) with eq. (1), (3) and (4) 
in Ref.^. 


B. Open system: 1 chemostat 

Introducing a chemostat breaks the concentration and 
mass conservation laws, but a new one arises (8). As 
a result, no affinity appears {s^ = 1, b = 1 and a = 
0) and the system evolves towards an equilibrium state 
compatible with the chemostat concentration Z'^y and the 
value of the component q (9) (the steady-state external 
current vanishes, I^y = 0). Also, since the system is now 
open, fcniax is infinite. 

Imposing the constraints on the expression of the steady 
state (18), namely 

< (1 - rf ( 20 ) 

[ Z'^^' = Z\r)'^y-\ 

we can express the variables Z^ and r numerically as 
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functions of q and , and obtain the equilibrium— 
exponential—distribution as a function of q and . 

C. Open system: 2 chemostats 

From two chemostats on, the infinite dimension of the 
system starts to play a role. As discussed in the previous 
section, two chemostats are not enough to drive the net¬ 
work towards a non-equilibrium steady state (s'*^ = 2, 
6=2 and a = 0): = 0 and = 0, where 

kyi and ky 2 represent the masses of the two chemostats 
(fcyi < ky 2 ). Thus, imposing the known values of the 
chemostat concentrations on the expression (18) leads to 

Z^yi _ 

( 21 ) 

Zky2 ^ ^i(y)fcy2-i^ 

which only admits physical solutions if > Z^^^. In 
this case, from (21) we obtain the equilibrium distribution 

fe-feyl 

^y2-^yl 

, ( 22 ) 

which is broader the smaller Z^^^ — is or the larger 
ky 2 — kyi is. When Z^^^ < Z^^^ the equilibrium concen¬ 
tration distribution becomes an increasing exponential 
which cannot be reached. As a result the system will 
enter a regime of continuous growth aimed at reaching 
that state (which we analyze in sec. V). 


where / G K determines the exact value of the fluxes. As 
for two chemostats, the set of equations in (23), in the 
variables Z^, r and J, does not exhibit physical solutions 
if the concentration of the largest chemostat is higher than 
the one of the smallest one, i.e. Z^^^ < Z^y^. On the 
other hand, whenever the above condition is not fulfilled, 
the stationary solution is unique and stable (appendix B). 
Solving the system (23) numerically, we obtain the values 
of Z^, r and I given Z^^b Z^y^ and Z^^^. In Fig. 4 
the distribution is shown for different values of these 
concentrations. 

The chemostat concentrations also determine the sign 
of the related fluxes: if the concentration of the second 
chemostat lies above the equilibrium distribution obtained 
by the first and third one, we have a continuous flow of 
mass from the intermediate chemostat towards the exter¬ 
nal ones (/ > 0, Fig. 4a). Vice versa, if the concentration 
of the second chemostat lies below the equilibrium distri¬ 
bution obtained by the first and the third one, we have 
a continuous flow of mass from the smallest and largest 
chemostats towards the intermediate one (/ < 0, Fig. 4b). 
Importantly, whatever physical value Z^yi, Z^^^ and Z^y^ 
assume, the system cannot exhibit a condition in which a 
net flux of matter from the largest species to the smallest 
one occurs. This is clear by looking at the jy^^-values in 
(14) used to express J^yy eq. (24): the sign of and 
are always the same, and opposite to the one of 


E. Open system: more chemostats 


Zk _ zkyi 


/ Z'"y2 \ 

I Z^ / 


D. Open system: 3 chemostats 

Three is the minimum number of chemostats able 
to drive the system in a non-equilibrium steady state 
(sec. IIIB). Indeed, a class of emergent cycles appears 
(s'^ = 3, 6 = 2 and a = 1) and the system exhibits a 
set of non-vanishing external currents. If we impose the 
values for the chemostats’ concentrations on the general 
expression for the steady state (18), we obtain: 


Going on adding chemostats, new independent classes 
of emergent cycles appear. The procedure for determin¬ 
ing the steady-state distribution is equivalent to that 
discussed is subsec. IV C and IV D. In these two cases 
we proved that when the largest chemostat has a con¬ 
centration greater or equal to that of the smallest one, 
the system does not reach a steady state. The same 
exact behavior has been observed numerically for more 
chemostats, hence we speculate that this property holds 
for an arbitrary number of chemostats. 


' ^fcyi 

Zky2 

Zkys 


Zi(y)feyi-1 

_ r^yl 1 (fi'\ky2 — kyl 

^l(y)fey2-l +-- 

K, 1 — r 

_ r^yl 1 — /'^'\^y3~^yl 

^i(y)fey3-i +- 

K, 

2. — ^^^^y3~^y2 

n 1 — f 


(23) 


1 — r 


As discussed in sec. IIIC, the external currents are 
subject to the same constraints as the emergent cycles, 
and can be written as linear combinations of them. Since 
we have one class of emergent cycles, characterized by the 
zy^y values in (14), we have that 


(24) 


As a final remark, we point out that the steady-state 
distributions do not depend on the value of the rate 
constant k. Indeed, solving the equations (20), (21), and 
(23) for Z^, f and / k, we obtain them as functions 
of the components and the chemostats’ concentrations. 
Since the latter do not depend on k, the same holds for 
Z^, r and / k. As a corollary is proportional to k 
and the same holds true for the entropy production (16). 


V. ASYMPTOTIC GROWTH REGIME 

We mentioned in the previous section that the system 
does not exhibit a steady state when the concentration 
of the largest chemostat exceeds that of the smallest 


/'=y' = 


1= 1,2,3 
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polymer mass, k 
(a) 



0 5 10 15 20 25 30 


polymer mass, k 
(b) 

FIG. 4: Non-equilibrium steady-state distributions for 
the system of monomer-exchanging polymers with 
three chemostatted species. In both of the plots, the 
chemostats—highlighted in green and by the arrows—are 
kyi = 2, ky 2 = 5 and kys = 10. The orientation of the ar¬ 
rows denotes the sign of the external fluxes of chemostats: 
arrows pointing up means chemostats leaving the system, 
i.e. > 0. The chosen chemostat’s concentrations are: 
plot (a) = 5, = 7 and = 2; plot (b) 

Z’^y^ = 5, Z^y‘^ = 1 and Z'^y^ = 2. 

one, Z^y^ < —we refer in the text to this con¬ 

figuration of chemostats leading to continuous growth 
as “unbalanced”. The dynamical fixed point moves out¬ 
side the region of physical solutions—namely to r > 1, 
see appendix B —and the system approaches the limit 
r —> 1. This indicates that the concentration of the single 
monomer species becomes negligible compared to the rest 
of the species. Hence the system grows towards an un¬ 
reachable steady state with an exponentially increasing 
concentration distribution. 

Fig. 5a shows the concentration distributions of an 
unbalanced system at different times before the numerical 
cut-off (more details are given in the related caption) is 
reached. These different distributions show that while 
the concentrations of the species between two chemostats 
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FIG. 5: Concentration distributions at different times are 
shown for system in unbalanced conditions. Different col¬ 
ors from red to violet correspond to exponentially increas¬ 
ing times. The set of plots is obtained by numerical solu¬ 
tion of the differential equation (5). Absorbing boundary 
conditions have been chosen, meaning that the concentra¬ 
tion at the cut-off—here set to fccut-ofi = 1000—is zero. 
We point out that this prescription is safe before the cut¬ 
off is reached. In plot (a) we report a system with three 
chemostats. The chemostat’s masses and the related con¬ 
centrations chosen are: Z® = 1, = 7 and Z®^ = 2. The 

concentrations of the species between the chemostats ba¬ 
sically overlap at times t > 1 and become steady. Beyond 
this time the growth only involves the species larger than 
the biggest chemostats. In plot (b) we consider a system 
with monomers and dimers chemostatted: Z^ = 3 and 
Z^ = 4. 


stabilize to steady values, the concentrations of the species 
larger than the biggest chemostat do not. Hence, the 
system continuously grows trying to populate the infinite 
size polymer. This behavior has been observed taking into 
account different number of chemostats and chemostat’ 
concentrations. 

In order to characterize this growth algebraically, we 
consider a system with monomer and dimer chemostats 































(fcyl = 1 and ky 2 = 2) such that (The typi¬ 

cal growth obtained numerically in this scenario is shown 
in Fig. 5b). Since the growth dynamics cannot be solved 
exactly, we assume that the asymptotic concentration dis¬ 
tribution can be parametrized by the (equilibrium) steady 
state expression (18) with time dependent parameters, i.e. 


for/fc>3. (25) 

where A(t) and a{t) are unknown real functions of time. 
To simplify the notation, let us denote the concentrations 
of the chemostats by = Z^^^ and = Z’^^^. The 
functions A{t) and a{t) can be determined by means of 
the differential equations for the total concentration Z 
and the total mass M: 


Z = = -kZ{Z^ - y2) _ 

M = -I^ - 2/2 = 

- kZ{2Z^ - 3y2 + yi) _ k2Y^Y^ + kY^Y^, 

(26) 

where, Z, M and the concentration of trimers Z^ assume 
the following form when the ansatz (25) is taken into 
account 


M{t) ~ + yi + 2Y^, 

(l-a(t)) 

Z^{t)^A{t). 


(27) 


When the equations are expressed in terms of A{t) 
and a{t), the stream plots for different values of the 
chemostats’ concentrations show that the ansatz captures 
the non-equilibrium phase transition occurring when the 
chemostats become unbalanced. Fig. 6. Indeed, for bal¬ 
anced chemostats, the system evolves towards a fixed 
point with a lying in ]0,1[, Fig. 6a. On the other hand, 
when the chemostats are unbalanced the fixed point lies 
beyond a = 1 signaling an asymptotic growth regime, see 
Fig. 6b. 

The numerical solution for A{t) and a(t) obtained using 
(26) and (27) accurately characterizes the asymptotic 
growth. Indeed, as seen in Fig. 7, when comparing the 
evolution of Z and M obtained from A{t) and a{t) with 
that obtained by solving numerically the rate equations, 
the former solution overlaps with the latter before the 
cut-off used in the numerics is reached. We find that the 
total concentration grows linearly with time whereas the 
mass quadratically. 

Taking into account the ansatz (25), the entropy pro¬ 
duction rate (16) becomes 


E ~ Rl^ In 


A{t) 

Y^{a{t)Y 


+ RI^ In 


A{t) 
Y^a{t) ’ 


(28) 


where, and can be written in terms of T^, Y^, 
A{t) and a{t) using eq. (26). The latter is plotted in 



(a) balanced chemostats 



0.0 0.2 0.4 0.6 0.8 1.0 

a(t) 

(b) unbalanced chemostats 


FIG. 6: Stream plot of the differential equations (26) ex¬ 
pressed in terms of the ansatz functions a{t) (abscissa) and 
A{t) (ordinate). When balanced chemostat concentrations 
are used, the fixed point lies for values of a{t) in ]0,1[: plot 
(a). The chemostats chosen for this plot are Y^ = 4 and 
y2 = 2. Vice versa, when the chemostats are unbalanced 
(y^ = 2 and Y^ = 4) the fixed point moves outside from 
the physical region (a(t) > 1): plot (b). 


Fig. 8, where it is compared with the numerical solutions 
for different cut-offs. The agreement with the numerical 
solution is not perfect but captures the linear asymptotic 
growth of the entropy production rate reasonably well. 
Also, we point out that the unbalanced dynamics shown 
in Fig. 8 exhibits an initial transient relaxation stage 
shown in inset. 

We conclude mentioning that the same ansatz could be 
used for systems characterized by more chemostats with 
unbalanced concentrations. Indeed, the growth always 
involves the species larger than the biggest chemostat, 
whereas the species between chemostats converge faster 
to proper steady values. Hence, fixing the concentration 
of these latter species, we could assume a growth like 
(25) for the species larger then the biggest chemostat and 
















9 



0 20 40 60 80 

time, t 
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FIG. 7: Total concentration (a) and total mass (b) as func¬ 
tions of time in the asymptotic growth regime. The nu¬ 
merical solution obtained using the ansatz (25) is plotted 
in green (dashed). These plots are compared with numer¬ 
ical solutions of the system of differential equations (5) 
with different cut-offs (blue curves). The chosen chemostat 
concentrations are: = 3 and = 4 while the initial 

condition imposed is Z^(t = 0) = |(|)*^. Finally, the cho¬ 
sen cut-off concentrations are: kc = 200 (dark blue curve), 
kc = 500 (blue curve) and kc = 1000 (light blue curve). 


perform the same analysis. 


VI. CONCLUSIONS 

This paper provides a kinetic description of systems 
made of glucans and processed by the class of enzymes 
known as D-enzymes. The action of the enzyme induces 
a monomer-exchange process between pairs of glucans 
which are distinguished by their mass or degree of poly¬ 
merization. Free monomers are not allowed to attach to 
other polymers^ implying that the total concentration 
and the total mass are conserved when the system is 
closed. The system’s dynamics is ruled by rate equations 
for the polymer concentrations endowed with mass action 
kinetics. We mimic physiological conditions by introduc¬ 



FIG. 8: Entropy production rate as a function of time in 
the asymptotic growth regime. The numerical solution ob¬ 
tained using the ansatz (25) is plotted in green (dashed). 
This plot is compared with numerical solutions of the sys¬ 
tem of differential equations (5) with different cut-offs 
(blue curves). In all the plot, the entropy production rate 
is given in units of R. The chosen chemostat concentra¬ 
tions are: = 3 and = 4 while the initial condition 

imposed is Z^(t = 0) = §(§)*’. The chosen cut-offs kc are: 
200 (dark blue curve), 500 (blue curve) and 1000 (light 
blue curve). Also, the inset shows in greater details the 
initial transient relaxation stage. 


ing chemostats which effectively describe the action of 
the environment by fixing the concentrations of certain 
glucans. In this scenario, chemostats represent species 
processed by the environment. For example, they may 
represent species which need to be processed and injected 
by the environment in the system, analogously, they may 
represent final products of the metabolic processes which 
are taken out of the system. Importantly, chemostatting 
the system amounts to open it and introduce driving 
forces on the non-chemostatted species. 

Our main results are summarized in Table I. We identi¬ 
fied three types of different long-time behaviors depending 
on the chemostatting conditions: equilibrium state, non¬ 
equilibrium steady state, and continuous growth of the 
system. The closed system as well as the open system 
with a single chemostat always relax to an equilibrium 
state. In presence of two chemostats the system will either 
relax to equilibrium or turn into a state of continuous 
growth depending on whether or not the concentration 
of the largest chemostat is lower than the concentration 
of the smallest one. We proved that this latter condition 
for growth holds true for up to three chemostats and 
conjectured that it is generally true based on numerical 
evidence. For more than two chemostats, if the concen¬ 
tration of the largest chemostat is lower than that of 
the smallest one, the system will reach a nonequilibrium 
steady state where the chemostats continuously exchanges 
matter across the system. Our results confirm that, even 
in the infinite-dimensional chemical network considered 
here, the number of chemostats equals to the number of 
broken conservation laws plus the number of emergent 















10 


cycles (see Tab. I). A proof of this equality for finite dimen¬ 
sional chemical networks is provided in Ref.^^. We also 
emphasized the role of the emergent cycles in driving the 
chemostatted chemical networks towards nonequilibrium 
steady states rather than equilibrium states^^. 

The metabolism of polysaccharides is a complex process 
involving many steps and several enzymes^ and its com¬ 
plete dynamical characterization is beyond the scope of 
the present paper. We focused on the dynamical character¬ 
ization of the disproportionating action of D-enzymes in 
the breakdown and synthesis processes of glucans^. Under 
physiological conditions, it has been pointed out that one 
of the possible role of D-enzymes in these processes is to 
produce glucans of large sizes (which are then processed 
by other enzymes) starting from medium sized ones^. 
Importantly, a production of glucose (monomers in our 
descriptions) is expected, too®. This disproportionating 
behavior can be reproduced in a (nonequilibrium) steady 
state by the three chemostats system depicted in Figure 
4a. The intermediate high concentration chemostatted 
glucans represent the species to be processed, while the 
low concentration chemostatted glucans represent the 
species to be produced—in this case the small and large 
glucans. In this scenario, a continuous flow of interme¬ 
diate glucans enters the system and consequently both 
the smaller and the larger glucans are steadily produced 
and expelled from the system (§ IVD). We stress that 
the production of the small glucans follows from the total 
concentration conservation law (§ III A), i.e. the fact that 
free monomers cannot attach to other glucans. As seen in 
sec. IV C, two chemostats are not sufficient to reproduce 
a nonequilibrium steady state. 

Also, under closed in vitro conditions, the equilibrium 
distribution (which has also been analyzed in Ref.® and 
can be equivalently obtained by means of Maximum En¬ 
tropy methods^) agrees with experiments^. This means 
that if chemostatting conditions could be implemented 
in vitro, our predictions could be verified experimentally. 
Such a procedure would also enable to engineer different 
polymer concentration distributions. 

The approach we developed could be easily extended 
to describe the behavior of more sophisticated forms of 
D-enzymes^ embedding further conservation laws. It is 
also relevant to study any type of exchange process or 
aggregation-fragmentation dynamics^® in an open system 
framework^®^^®, emphasizing the importance of conserva¬ 
tion laws and providing more insights into the mechanisms 
driving these processes out of equilibrium. 
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Appendix A: Steady-state distributions 


The generic expression for the steady-state distribu¬ 
tion (18) can be obtained as follows. The steady-state 
equations can be expressed as 

Z {Z'^ - Z^} + = 0, 

Z - 2Z'= -b Z^-^} + Z^ {Z'= - Z'^-^} 

jk 

= —4fc„GnY, for fc > 2. 

Hi ^ 

(Al) 

Defining the variable AZ^ = Z^ — Z^~'^ they become 

ZAZ'^ + Z^Z^ = 0, 

Z {AZ’^+^ - AZ’^} + Z^AZ^ |-^2) 

jfc 

= —for A: > 2. 

Hence, by hierarchically substituting these expression one 
into the other, and using the variable f =1 — Z^/Z, we 
obtain 


AZ'^ = -(1-f) Z^f 

Py 




-yfc-V-i0(fc_A:y_i) ^ (A3) 


ky G Oy 


where ©(•) represents the discrete step function: 


Q{k) = 
Finally, 


0 if fc < 0 , 
1 if fc > 0 . 


(A4) 


Z'^ = ^ AZ* = Z^{r)'^-^- 


2=1 


+ E — ep-ky-i), (A5) 

fcy G^Y 

which corresponds to the equation (18) in the main text. 


Appendix B: Three chemostats steady state 

We discuss the uniqueness and stability conditions for 
the steady state when three chemostats are present. 

From the constraints on the steady state (23) and from 
the condition for the external currents (24), we can write 
a single steady state condition involving just f as variable: 

_|_ pyiz'^y^'^ (^^i^Gi+,y'‘y 3 _ 

(jyfcyl +iyU3)^fey2(^)Uyl_ 

_ (Bl) 

('j,fcy3^fcyl + ,yUl^fcy3) (b)'^'‘^"-b 
(jyfcy3^fcy2 _p,yfeyl^fcy3) = Q. 
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• Xq < 1 and 2/0 < 1 (^o < 0). This condition im¬ 
plies the following configuration for the chemostats: 

In this case we have always one and only one solu¬ 
tion. Indeed, the center of the hyperbola (xo,yo) 
takes (0,1) x (0,1), and the upper right branch of 
the hyperbola always intersects the power law in 
X = 1 (which is non-physical). The left lower one, 
instead, always intersects the power law for values 
in (0,1) X (0,1) since yh(x = 0) > 0 (fig. 9a). 

• xo < 1 and_ 2 /o > l_(^o > 0)._This condition corre¬ 
sponds to > Z^y^ and Z^y^ > Z^y^. 

In this case we have one solution if and only if 
Z^y^ > Z^y^. The center of the hyperbola lies in 
(0,1) X (1, oo) and the upper left branch of the hy¬ 
perbola never intersects the power law. The right 
lower one, instead, always intersects the power law 
in a: = 1, y = I (fig. 9b). We have a further 
intersection in the physical region if and only if 
> ^ , which holds iff Z^y^ > Z’^y^— 

X — 1 X — 1 

indeed, x* : yh{x*) = 0 is such that x* > 0, for any 
choice of the chemostats. 


dyp 

dtc 


FIG. 9: Plots of the hyperbola (dark purple curve) and 
power law (light purple curve) in (B2) for different con¬ 
figurations of parameters. The center of the hyperbola is 
highlighted by a dark purple dot, while the physical region 
by the dashed orange lines. 


Let us define the variables x = (r)'' and y = (f)’' , 

so that the above-expressed steady-state condition can be 
written as the intersection of two curves: a rectangular 
hyperbola and a power law function 



(B2) 


where the coefficients are given by 

(j/^yl _|_ jy^ys) Z^y^ 

Xo = 


yo = 


zo = 


l/^y3 Z^y^ t/^yi ^^y2 

jyfcyS^fcyl _|_ jyfeyl ^^yS 


l/^y3 Z^y^ i/^y^ Z^y^ 

iykyljyky3 (^2^y2 _ ^^yl ) (^Z'^y^ - ) 

^l/ky3 Z^yl _|_ l/kyl ^ky2'^ 


(B3) 


[The subscripts h and p simply help us to distinguish the 
two functions.] From a geometrical point of view, physical 
solutions are represented by those intersection points lying 
in {x,y) G (0,1) x_(0,1)._In order to prove that this 
happens whenever Z^y^ > Z^y^ we observe that all of the 
possible configurations of chemostat concentrations are 
described by the following four cases for the parameters 
xo and yo- 


• a;o > I and yo < I [zq > 0). This condition corre¬ 
sponds to: Z'^y^ < Z^y^ and Z'^y^ < Z^y^. 

Once again, we have one solution if and only if 
Z^yi > Z^y^. The center of the hyperbola lies 
in (I,oo) X (0,1) and the right lower branch of 
the hyperbola never intersects the power law. The 
upper left one, instead, always intersects the power 
law in a; = I, y = I (fig. 9c). We have a further 
intersection in the physical region if and only if 
> ^ , which holds iff Z^y^ > Z’^y^— 

X—1 X—1 

indeed, yh(0) > 0 for any choice of the chemostats. 

• a;o > I and yo > I (^o < 0). This condition im¬ 
plies the following configuration for the chemostats: 
Z^yi < Z^y'^ < Z^y^. 

In this case we have no solutions. Indeed, the center 
of the hyperbola lies in (a;, y) G (1, oo) x (1, oo) and 
neither the upper right nor the lower left branch of 
the hyperbola intersects the power law in the physi¬ 
cal region. The left lower one, indeed, always inter¬ 
sects the power law in (1,1) which is non-physical 
(fig. 9c). 

Summarizing, we have a unique steady state whenever 
the concentration of the largest chemostat is higher than 
the concentration of the smallest one: Z^y^ > Z^y^. 

Stability. In order to prove the stability of the fixed 
point we resort to the following Lyapunov function: 

L = ^Z'=ln |J-(Z-Z,). (B4) 

k ® 

It is easy to prove that this function is always positive 
and vanishes only for = Z^, where Z^ represents the 


dyp 

dx 
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steady-state solution. If the steady-state solution exists, 
namely if exists ■. = 0, the time derivative of the 

Lyapunov function (B4) can be written as 


dt 





(B5) 


Close to the steady state the above derivative is negative. 
Spanning the phase space with small perturbations on 
every concentration, we always obtain ^ < 0, where the 
equal sign is reached only at the steady state. Disregard¬ 
ing the infinite dimension of the phase space, we consider 
the independent set of perturbations labeled with the 
index k'^ and quantified by the small real value e 

^ g^fcifex^ |g| ^ jnin(B6) 


Embedding these perturbation in (B5) and using the rate 
equations (5) we obtain 


dL 


dt 


dL 


dt 

ryk'^ 

•^S 






2Z^ + 2Z 


-z^- ^-Z} 


for = 1, 


(B7) 


) for k'^ yf 1, 


which are always negative, no matter the sign of the 
perturbation. 
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